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We explore the sensitivity of limb darkening coefficients computed from stellar 
atmosphere models to different least-squares fitting methods. We demonstrate 



> 



O 
Oh 



CO 

c3 



that conventional methods are strongly biased to fitting the stellar limb. Our 



suggested method of fitting by minimizing the radially integrated squared resid- 



ual yields improved fits with better flux conservation. The differences of the 
obtained coefficients from commonly used values are observationally significant. 
\q | We show that the new values are in better agreement with solar limb darken- 

ing measurements as well as with coefficients reported from analyses of eclipsing 
binary light curves. 
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"*-! ■ 1. INTRODUCTION 



With the exception of our Sun, up to recently the observational significance of stellar 
limb darkening had been limited to the analysis of eclipsing binary light curves (Popper 1984). 
A range of new observational techniques and strategies sensitive to limb darkening have been 
developed since then. These include for example stellar surface interferometry (Quirrenbach 
et al. 1996), observation of gravitational microlensing caustic-crossing events (Witt 1995) 
and extrasolar planet transits (Deeg, Garrido, & Claret 2001). When analyzing any such 
measurements, expected intensity profiles from stellar atmosphere models are often used in 
the process of determining other event parameters. However, with sufficiently good quality 
data these observations can also be used for direct measurement of stellar limb darkening. 



Such measurements are usually done by assuming simple analytical limb darkening laws. 
The measured coefficients of these laws are then compared with their theoretically predicted 
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values, which are obtained by fitting the laws to model atmosphere intensity profiles. How- 
ever, these predicted values depend significantly on the adopted method of fitting (e.g., 
Diaz-Cordoves, Claret, & Gimenez 1995). Here we demonstrate that even variants of the 
commonly used least-squares methods can yield very different coefficient values, and suggest 
a preferred method of fitting. 

In the following §2 we discuss various fitting methods and introduce our suggested 
method. We study the fit qualities of selected methods in §3 and demonstrate the differences 
between obtained coefficients on the example of ATLAS model atmospheres (Kurucz 1993a). 
In §4.1 we compare these coefficients with measured solar limb darkening, in §4.2 with limb 
darkening coefficient estimates available from eclipsing binary light curve analyses. We 
discuss variations of the suggested fitting method in §5 and summarize the main conclusions 
in §6. 



2. LIMB DARKENING LAW FITTING METHODS 

One of the main results obtained from stellar atmosphere models is the center-to-limb 
variation of specific intensity I v of the emitted radiation field computed for a range of wave- 
lengths. The center-to-limb dependence is obtained by computing radiative transfer along a 
set of rays emerging from the atmosphere under different viewing angles 6. Broadband limb 
darkening of the intensity / can then be computed from l v for different photometric filter 
systems. Intensity values are traditionally parametrized by /i = cos#, ranging from 1 at 
stellar disk center to at limb, or alternatively by the radial position on the disk r = sin 9, 
ranging from at center to 1 at limb. 

Various analytical laws have been suggested for approximating stellar limb darkening, 
most of them in the form of linear combinations of simple functions of fi. The most basic 
limb-darkening law, generally referred to as "linear limb darkening" , is given by 

I L {y) = I [l - u(l - //)] , (1) 

where the linear limb-darkening coefficient u determines the shape of the limb-darkening 
profile, and Iq is the intensity of the linear law at disk center. Other laws have been intro- 
duced to get a more accurate analytical description, by adding a third term to equation (1), 
for example the quadratic, square-root, or logarithmic laws (see Claret 2000). Claret (2000) 
also proposed a more complex law which is a linear combination of five terms. 

A range of methods have been used for the seemingly straightforward task of computing 
limb-darkening coefficients from the intensity profiles of model atmospheres, each of them 
producing generally different values. 
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2.1. 



PREVIOUS METHODS 



The two main groups of fitting methods are least-squares based methods (e.g., Manduca, 
Bell, & Gustafsson f 977; Claret & Gimenez 1990; Diaz-Cordoves et al. 1995; Claret, Diaz- 
Cordoves, & Gimenez 1995), and methods requiring flux conservation (e.g., Wade & Rucinski 
1985; Van Hamme 1993). The primary motivation for the latter methods is the observation 
that by concentrating on the quality of the fit of the profile shape one often ends up with 
a low accuracy of flux conservation. This problem is pronounced especially in the case of 
the linear law, less so in higher order laws. Computation of higher order coefficients in the 
flux conservation methods requires additional constraints, for example conservation of mean 
intensity (Van Hamme 1993), or fixing the intensity value at fx — 0.1 close to the limb (Wade 
& Rucinski 1985). 

Least-squares based methods compute limb-darkening coefficients mostly by minimizing 



where is the intensity of the model atmosphere at disk position /ij, and II is the limb- 
darkening law of interest. In support of the least-squares methods Claret (2000) argues that 
flux-conservation methods fail in the main reason for using analytical laws, which is to get a 
reasonable fit of the shape of the intensity profile. A good fit to the profile - particularly in the 
case of higher order laws - should inherently lead to sufficiently accurate flux conservation. 
In comparison, by concentrating on integrated quantities flux-conservation methods largely 
ignore the shape of the profile. 

We note here that in earlier works based on either of the methods (e.g., all those cited 
in the first paragraph) the central intensity Iq was kept fixed at the model atmosphere value 
and only the shape-defining limb darkening coefficients were fitted. However, there is no a 
priori reason to fix one particular point in the intensity profile at its exact value and fit only 
the remaining points. On the other hand, this approach leads to an overall worse fit. In 
more recent works this constraint has been apparently dropped (Claret 2000; Barban et al. 
2003). In the case of the linear law from equation (1) the simplest approach is to fit the 
linear parameters I and ul and compute the limb darkening coefficient u as the ratio of 
the obtained values. All the other laws mentioned above can be fitted in a similar manner. 

The least-squares method has its own further ambiguities. As pointed out by Diaz- 
Cordoves et al. (1995), the obtained coefficients depend on the distribution of the rays /ij 
provided by the model atmospheres. For example, in the commonly used Kurucz ATLAS9 
models (Kurucz 1993a,b,c, 1994) the spacing of the 17 computed rays decreases sharply 
at low values of //, close to the limb. A straightforward fit minimizing equation (2) then 
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gives a high weight to the quality of the fit close to the limb on account of the fit quality 
closer to the center. In order to get a less biased fit Dfaz-Cordoves et al. (1995) as well as 
subsequently Claret et al. (1995) or Claret (2000) used only a subset of the available rays 
more evenly spaced in fi. Similarly, Barban et al. (2003) who used their own ATLAS-based 
model atmospheres computed 20 rays evenly distributed in /i. 

2.2. SUGGESTED METHOD 

A simple way to overcome problems caused by the spacing of discrete points \pi, Ij\ is to 
suitably interpolate the intensity between the points, so that we get a continuous function 
which passes through all the model points, /(/ij) = For this purpose we can use 
cubic spline interpolation. The Jj versus fii dependence is typically smooth, without large 
jumps or wiggles, so that such a spline gives a reasonable interpolation. The limb darkening 
coefficients can now be computed by minimizing the integrated squared residual 

K = j[h-I] 2 dfi, (3) 

where all integrals appearing in the least-squares solution can be evaluated to arbitrary 
precision. The result obtained in this way is what previous authors were trying to achieve by 
selecting regularly spaced points. In the suggested approach none of the points obtained from 
the model atmosphere computations are omitted, thus no useful information is neglected. 
The spline interpolation provides a sensible guess about the behavior of the intensity in 
between these points. We discuss the sensitivity of the results to the particular choice of 
spline interpolation in §5. 

However, the ambiguity due to the spacing of points is overcome only seemingly. The 
new ambiguity lies in the choice of the integration variable in equation (3). In this example, 
integration over /i gives equal weight to the squared residual in equal-sized intervals of 
fi. For instance, the squared residual in the radial interval from the center to r = 0.87 
(corresponding to /i — 0.5) has the same weight as in the remaining interval close to the 
limb. By transforming from jj, to the radial coordinate r in the integral in equation (3) we 
get 

Aj = J[I L -I] 2 ril-r 2 y l ' 2 dr. (4) 

We see that the residual at the disk center has zero weight and the weight function grows 
monotonously towards the limb where it diverges. The choice of fi as the integration variable 
thus still gives results strongly biased to fitting the limb of the star. In order to get a 
reasonable fit even closer to the center, we decide to integrate over r and compute limb 
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darkening coefficients by minimizing 



A 2 




(5) 



The suitability of the choice of integration variable (or equivalently weighting) should 
be primarily demonstrated in confrontation with observational data, and we do so in §4. 
In support of our choice we briefly consider here some relevant astrophysical situations, 
namely, eclipsing binaries and caustic-crossing microlensing events. In either of these cases 
the weighting kernel that is convolved with the intensity to obtain the flux is time-dependent, 
depending on the phase of the binary or the position of the passing lens, respectively. Hence 
neither scenario directly provides a single ideal weighting. However, if we look for simplicity 
at the case of a 90° inclination binary, the radial position of the occulting component shifts 
linearly with time. Moreover, for example in the case of annular eclipses the occulted area is 
of constant radial size. Similarly, in zero impact parameter single-lens microlensing the radial 
position of the lens changes linearly with time, at each position giving maximum weight to 
the point directly behind the lens. Such considerations suggest that equal radial weighting 
is a natural and reasonable choice. 

In the following section we illustrate the differences between the results obtained by the 
discussed fitting methods on the example of widely used model atmospheres. 



We demonstrate the different fitting methods on the ATLAS9 model atmospheres com- 
puted by Kurucz (1993a,b,c, 1994) using the intensity files available on Robert Kurucz's 
website 1 . The models in the parameter grid have metallicities ranging from [Fe/H] = +1 to 
-5, effective temperatures from T e ff = 3500 to 50000 K, surface gravities in CGS units from 
log g = to 5, and microturbulent velocities for solar metallicity models from Vt — to 8 
kms" 1 , for other metallicities only 2 kms^ 1 . For each model specific intensity I u is given 
for wavelengths from A = 9.09 nm to 160 /im spaced by 2nm in the optical range, and for 
17 rays passing through the atmosphere: 10 at /i values from 1 to 0.1 spaced by 0.1, and 7 
additional rays close to the limb at 0.25, 0.15, 0.125, 0.075, 0.05, 0.025, 0.01 (see Figure 1). 
We obtain intensity profiles for specific photometric filters by integrating 
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DEMONSTRATION ON ATLAS ATMOSPHERES 
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1 http://kurucz. harvard.edu 
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where c/A 2 is the specific intensity conversion factor (from per frequency to per wavelength), 
and S(X) is the filter response function. In this work we use the standard BVRI filter 
response functions given by Bessell (1990). We formally extend the intensity profile to 
fi = by linear extrapolation from the two points closest to the limb, and obtain /(//) by 
interpolating the profile using cubic splines with natural boundary conditions. We discuss 
the influence of different boundary conditions in §5. 

We compare four different least-squares fitting methods. First, we fit all 17 points 
provided by Kurucz's models minimizing \ 2 given by equation (2). Second, we use the 
same equation to fit a subset of 11 points with 10 regularly spaced /i values from 1 to 0.1 
plus \i = 0.05 - these are the points fitted by Dfaz-Cordoves et al. (1995), Claret et al. 
(1995) or Claret (2000). Third, we fit the spline-interpolated intensity / minimizing A 2 
from equation (3) using integration over /i. Fourth is our suggested method of fitting / by 
minimizing A 2 from equation (5) using integration over r. 

In Figure 1 we present as an example the V band profile of a T e ff = 4250 K, log g = 4.5, 
[Fe/H] = —3.5 model atmosphere plotted separately as a function of r and /i. The plots 
include the original points and their spline-interpolated intensity profile / (here normalized 
to unit r— integrated rms intensity), as well as linear limb darkening laws fitted by the four 
methods and their residuals from /. As expected and as pointed out by Diaz-Cordoves et 
al. (1995), the 17-point fit has the heaviest bias to the limb. The 11-point fit here closely 
approximates the mu— integrated fit which has the lowest m-u-integrated squared residual 
(see lower right panel), and both are only slightly less biased to the limb. The r— integrated 
fit by definition has the lowest r— integrated squared residual (see lower left panel). The 
non-integrated residuals of this fit are overall the smallest from the center to roughly jj, = 0.2 
(i.e., r = 0.98). Although the quality of the //—integrated fit is substantially better on 
the remaining interval, the radial width of this interval is only 0.02. We conclude that the 
r— integrated fit gives the best linear limb darkening approximation of the intensity profile. 

In order to compare the fitting methods more generally, we fitted BVRI profiles of 
the full range of models described above, a set comprising altogether 38324 profiles of 9581 
model atmospheres. Besides the linear law we used quadratic, square- root, logarithmic, and 
Claret limb-darkening laws (for definitions see Claret 2000). For each fit we computed two 
parameters as measures of fit quality. The first is the relative rms residual defined as the 
ratio of rms residual to rms intensity 

A r 



(7) 



-7- 

Our second parameter is the relative flux excess of the obtained fit over the original profile 

AF _ J I L r dr 
F J lr dr 

The average and maximum absolute values of these parameters for the five limb darkening 
laws and four fitting methods are presented in Table 1. Judging by the residuals a, for all the 
laws r— integration by definition gives the best results (for the linear law average o = 1.25%, 
worst fitted profile has a = 3.53%), followed by the substantially worse 11-point fit (2.34%, 
5.55%), /i-integration (2.71%, 6.28%), and the poorest 17-point fit (3.33%, 6.99%). 

Comparing the fitting methods by relative flux excess yields a different picture. For 
all the laws \x— integration gives by far the best flux conservation (here limited by machine 
precision, for the linear law formally average \AF/F\ = 1.18 x 10~ 16 , the highest value for 
a profile |AF/F| = 5.96 x 1(T 16 ), followed by r-integration (2.27 x 1(T 3 , 8.65 x 1(T 3 ). 
Except in the case of the linear law where the order of the remaining methods is reversed, 
these are followed by the 11-point fit (4.38 x 1(T 3 , 1.29 x 1(T 2 ) and the 17-point fit (3.63 x 
10~ 3 , 1.42 x 10~ 2 ). The reason for the excellent flux conservation 2 in fj,— integration is the 
similar weighting of the intensity in the flux integral, as seen in equation (8). The factor r 
gives zero weight at the disk center and highest (though not divergent) weight at the limb. 
Nevertheless, in view of the roughly two times higher average residual in comparison with 
r— integration for the four simpler limb darkening laws, we conclude that this method suffers 
from a similar problem as the non-least squares flux-conservation methods (Wade & Rucinski 
1985; Van Hamme 1993), i.e., achieving good flux conservation at the cost of a poor fit of the 
intensity profile shape. The obtained results support the use of r— integration fitting of the 
spline-interpolated intensity profile as the method of choice for computing limb darkening 
coefficients. 

We note that Table 1 can be also used for example to compare the overall quality of the 
limb darkening laws for the different fitting methods. However, discussions of the different 
analytical laws are beyond the scope of this paper. 

The main result of this section for practical purposes is the numerical difference between 
the limb darkening coefficients produced by different fitting methods. We limit ourselves here 
to comparing those obtained by our suggested method of r— integration (hereafter u r ) with 
those obtained by the 11-point fit (hereafter %), the most frequently used coefficients (e.g., 
Diaz-Cordoves et al. 1995; Claret et al. 1995; Claret 2000). For illustration we plot BVRI 



2 Note that accuracy better than ~ 10 4 exceeds the flux accuracy of the underlying model atmosphere 
data, as shown in §5. 
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values of u r in Figure 2 as a function of T e ff. In order to avoid crowding we include only solar 
abundance models with v t = 2kms~ 1 . The specific range of log g of the models depends on 
T e ff. The main tendencies are well known: overall the coefficient decreases with increasing 
temperature and increasing wavelength. The rms residuals and flux excesses for these models 
and both types of fitting are shown in Figure 3 as a function of T e ff. The r— integrated fits 
have lower residuals and better flux conservation, in agreement with the results in Table 1. 

Figure 4 illustrates the difference between u r and Uu for the full range of model atmo- 
spheres, plotted as a function of The coefficient difference is substantial, ranging from 
-0.14 to 0.11 with most of the u r lower than Uu, corresponding to flatter limb darkening 
profiles. At the same time, the most limb darkened profiles are even more limb darkened 
when using u r . The results are plotted again in Figure 5 as a function of T e ff for each 
photometric band separately. The greatest variation occurs for lower temperature models 
(T e ff < 7000 K) and is much reduced for higher temperature models (T e ff > 10000 K). The 
range of coefficient differences decreases from the B band to the / band. We note specifi- 
cally that 11-point fitting overestimates the limb darkening coefficient in the I band for all 
computed models, in the R band for all models with T e ff > 4500 K, in the V band for all 
models with T e jj > 6000 K, and in the B band for all models with T e ff > 7000 K. At lower 
temperatures coefficients are overestimated for some profiles and underestimated for others. 

A comparison of individual coefficients of higher order limb darkening laws reveals there 
can be even greater differences between the values obtained by different fitting methods. 
Nevertheless, as seen from Table 1, at the same time the overall fit quality is better, thus 
the obtained combinations of coefficients produce more similar limb darkening shapes. 

4. COMPARISON WITH OBSERVATIONAL RESULTS 

Following the discussions of relative merits of different limb darkening law fitting meth- 
ods we turn to the question of observational relevance. We note that the differences between 
the linear coefficients demonstrated in the previous section are large enough to be detectable 
not only in solar limb darkening measurements, but for example also in observations of 
caustic-crossing microlensing events or eclipsing binaries. In Cassan et al. (2006) we showed 
that the coefficients measured in caustic-crossing microlensing events involving K giants 
in the Galactic Bulge are in better agreement with the values obtained by r— integration 
than with the 11-point fit values of Claret (2000). Here we present a comparison of these 
coefficients using solar limb darkening measurements and eclipsing binary results. 
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4.1. SOLAR LIMB DARKENING 

For our comparison we use measured solar limb darkening data from Cox (2000, p. 356) 
for the four best-sampled wavelengths in the optical range: 450, 500, 550, and 600 nm. We 
compare these limb-darkening profiles with linear limb darkening fits for the same wave- 
lengths of the Kurucz (1994) solar model with T eff = 5777 ^, log g = 4.4377, [Fe/H] = 0.0, 
v t = 1.5 kms~ l . The r— integrated and 11-point fits give us coefficients u r and uu, respec- 
tively. We take these fixed coefficients and fit only the central intensity I from equation (1) 
to get a best fit to the spline-interpolated solar limb darkening. 

The relative residuals of these fits are plotted in Figure 6 for each of the wavelengths 
as a function of r. The values of the relative rms residual a and relative flux excess AF/F 
are given in each of the panels. We see that for all four wavelengths u r gives better results 
than Mn, its a is lower by a factor of 1.3 - 1.7 and AF/F is lower by a factor of 5 - 21. In 
all cases "Un yields a much steeper limb darkening profile than u T and the measured profile. 
We conclude that the theoretical limb darkening coefficient obtained by r— integrated fitting 
of the model atmosphere data is in substantially better agreement with the measured solar 
limb darkening than the coefficient obtained by 11-point fitting. 

4.2. LIMB DARKENING FROM ECLIPSING BINARY ANALYSES 

Obtaining limb darkening coefficients from the analysis of eclipsing binary light curves 
is no simple task. Already Popper (1984) estimated that getting a limb darkening coefficient 
with a precision of 0.1 or better requires at least 100 light curve points within minima 
with a 0.005 mag photometric accuracy. Not surprisingly, there are only a few eclipsing 
binary systems for which limb darkening coefficients have been extracted and published. 
These limb darkening measurements have been obtained by different teams using different 
methods, approximations, and approaches - and are therefore also weighted by different 
problems. In this section we do not re-analyze any of the binary light curves. We merely 
perform the simplest test for our purposes, collating the published values and checking their 
agreement with theoretical coefficients obtained from model atmospheres as described above. 

More specifically, we searched the literature for eclipsing binaries with measured limb 
darkening coefficients, for which estimates of T e ff, logg, and [Fe/H] were also available 
(we note here that these estimates may sometimes also be problematic). In a number of 
these binaries the coefficients for the primary and secondary stars were not measured inde- 
pendently. Instead, their difference was kept fixed at a value expected from some specific 
model atmosphere coefficients. We restrict our sample to the systems analyzed without this 



-10- 



constraint plus to those for which the fixed difference values agree within 0.01 with those 
expected from our u r coefficients as well as Claret (2000) coefficients (hereafter uq)- The 
measured BVR linear limb darkening coefficients for these nine binaries are listed in Table 2 
together with the corresponding values of uc and u r . These were obtained by trilinear in- 
terpolation of the values from Claret (2000), and of the values computed for the Kurucz 
ATLAS9 model atmosphere grid (Kurucz 1993a,b,c) by r— integration fitting, respectively, 
using the observed stars' T e ff, log g, and [Fe/H] values. 

For five of the binaries (BS Dra, WW Cam, V459 Cas, WW Aur, RW Lac) all measured 
coefficients are closer to the u r values, for MU Cas they are closer to the uc values, while 
for the remaining three (EE Peg, FS Mon, GG Ori) some are closer to u r and others to uq- 
Looking at the latter four systems, in the case of MU Cas we note that Lacy, Claret, & 
Sabby (2004b) first obtained u = 0.32 for both components in an unconstrained fit - a value 
in perfect agreement with u r . The value quoted in the table was obtained by fitting keeping 
a fixed luminosity ratio. Turning to EE Peg (Lacy & Popper 1984) we see that u r values for 
both components are also in agreement with the measurements. The FS Mon coefficients 
(Lacy et al. 2000) have no error bars 3 and the agreement with u r in the V band is only 
marginally worse than with uc, while the agreement with uc in the B band is substantially 
worse than with u r . In the last system, GG Ori (Torres et al. 2000), the B and R values agree 
better with u r while the V coefficients agree better with uc- Here the situation is similar 
to MU Cas: the original unconstrained fits from two datasets produced V band coefficients 
0.48±0.08 and 0.41±0.05, values better consistent with u r = 0.45 than with u c = 0.51. The 
unconstrained R band value of 0.33 ± 0.13 is also in much better agreement with u r = 0.38 
than the value quoted in the table. 

The limb darkening of eclipsing binaries in which the radii of the components are a 
sufficiently large fraction of the orbital radius (> 10 - 15%) may be influenced by the reflection 
effect. In such systems the measured limb darkening coefficient may not agree with the value 
given by the star's T e ff. In order to avoid this potential effect, we check the results for binaries 
in our selection with the smallest radii. These are RW Lac, V459 Cas, and GG Ori, with 
radii in the range 4 - 7%. From the values in Table 2 and the discussion of GG Ori above 
we see that the conclusions about the comparison of the two coefficients remain valid. 

Just as in the case of caustic-crossing microlensing events, our conclusion is that for 
eclipsing binaries limb darkening coefficients obtained by our suggested method give over- 
all better results than the widely used Claret (2000) coefficients. The systems in which 
original unconstrained fits give better agreement with expected values than photometrically 



3 This is also the case with RW Lac. 
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adjusted fits indicate possible problems caused by poor flux conservation or more generally 
the inadequacy of the linear limb darkening law. 

5. DISCUSSION 

When comparing the fitting methods in §3 one could also evaluate the relative \i— integrated 
rms residual, differing from a in equation (7) by using A M in the numerator and integrating 
over \i instead of r in the denominator. Such a residual is by definition lowest for the results 
obtained by fi— integration. For the linear law we get an average residual 3.28% and max- 
imum residual 6.94% - both values higher than the corresponding r— integrated residuals. 
We get similar results even for the second order laws. This only demonstrates that fitting 
of simple analytical laws is more difficult in terms of fi than r. As noted before, fitting by 
integration over fi gives very high weight to the limb, where attempts to fit the more complex 
intensity dependence go on account of the quality of the remaining fit. Only Claret's 5-term 
law has an average //—integrated residual lower than its average r— integrated residual, even 
though its maximum /i— integrated residual remains higher. This law therefore tends to fit 
the limb darkening close to the limb slightly better than closer to the center. 

In this work we interpolated intensity data points using cubic spline interpolation in 
terms of \i with natural boundary conditions, i.e., zero second derivative at \i = and /i — 1. 
In order to check the sensitivity of the results to the interpolation method, we tested using 
cubic spline interpolation in terms of r with zero first derivative at center and zero second 
derivative at limb. We note that this interpolation is not esthetically optimal, as in the 
region closest to the limb the r dependence is very steep, which can produce small wiggles 
in the spline between the last points, noticeable when plotted as a function of fi. Still, the 
obtained linear coefficients differ at most in the third decimal place. 

Similarly we tested the influence of different boundary conditions for the fi spline. 
Changing the boundary condition at the limb or the method of extrapolating the limb inten- 
sity had a negligible influence, on the order of 10~ 7 in flux or less. Changing the boundary 
condition at the center to equality of first and second derivatives (which gives an interpo- 
lation equally pleasing to the eye) had a larger influence, though still only on the order of 
1CT 4 in flux or less. The difference is due to the larger spacing of model atmosphere rays at 
the center. Based on this test we can consider the value of 1CT 4 as a limiting flux accuracy 
for the model atmosphere data used in this work. 
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6. SUMMARY AND CONCLUSIONS 

We have studied the sensitivity of limb darkening coefficients computed from stellar 
atmosphere models to different least-squares fitting methods. Tests with Kurucz ATLAS9 
model atmospheres (Kurucz 1993a,b,c, 1994) show that best results are obtained by suitably 
interpolating the available intensity points and fitting limb-darkening laws to the obtained in- 
tensity profile I by minimizing the r— integrated squared residual introduced in equation (5). 
Other methods give too much weight to the quality of the fit in a narrow region at the stellar 
limb. 

The differences between the coefficient values computed by the different methods are 
substantial. In particular, we have shown that linear limb darkening coefficients of BVRI 
profiles of Kurucz models computed by the suggested method differ by as much as 0.14 from 
values obtained by a commonly used fitting method with points roughly evenly spaced in 
terms of /j,. Such differences are observationally significant. The plots in Figure 6 illustrate 
that coefficients computed from a solar model atmosphere by the preferred fitting are in 
better agreement with solar limb darkening data. In Table 2 we have further demonstrated 
that available limb darkening measurements from eclipsing binaries are more compatible 
with the newly computed coefficients than with those from Claret (2000). We have reached 
a similar conclusion previously with limb darkening measurements from caustic-crossing 
gravitational microlensing events (Cassan et al. 2006). 

Other methods of measuring limb darkening such as stellar interferometry or obser- 
vation of extrasolar planet transits can use a similar approach when comparing measured 
coefficients with theoretical predictions. However, all these methods may run into problems 
with analytical limb darkening laws. The limb darkening signature in the data is usually not 
strong enough for a reliable extraction of more than a single limb darkening coefficient, i.e., 
the linear coefficient. At the same time, even at the photometric accuracy of the available 
data the linear law is in many cases an inadequate approximation, largely due to its poor 
flux conservation and generally poor matching of the limb darkening shape. 

The drawbacks of measuring coefficients of analytical laws can be bypassed by using an 
alternative description of limb darkening, obtained by principal component analysis of an 
ensemble of model atmosphere intensity profiles (Heyrovsky 2003). At the cost of not being 
analytical, such an approach gives the best possible limb darkening description as a linear 
combination of the least number of terms. The principal component description is therefore 
also ideally suited for measuring stellar limb darkening, as demonstrated by Heyrovsky (2003) 
in the case of microlensing events. 
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Table 1. Limb darkening fit quality for different least-squares fitting methods 



LU law 


Met nod 


Average a 


Max. a 


Average \l\r / r 


A/Tr,->r 1 A TP 1 T?\ 

Max. | I\r / r \ 






{%) 


{%) 






Linear 


17-point 


3.327 


6.987 


3.63 x 10~ 3 


1.42 x 10~ 2 




11-point 


2.344 


5.554 


4.38 x 10~ 3 


1.29 x 10~ 2 




/x— int. spline 


2.706 


6.275 


1.18 x 10~ 16 


5.96 x 10~ 16 




r— int. spline 


1.253 


3.527 


2.27 x 10- 3 


8.65 x lO" 3 


Quadratic 


17-point 


1.180 


3.303 


3.85 x 10~ 3 


1.07 x 10" 2 




1 1 -point 


0.653 


1.785 


1.10 x 10~ 3 


2.98 x 10~ 3 




fx— int. spline 


0.980 


2.609 


6.38 x 10~ 16 


1.75 x 10~ 15 




r — int. spline 


0.417 


0.991 


2.35 x 10 -4 


6.15 x 10~ 4 


Square-root 


17-point 


0.498 


3.545 


1.13 x 10~ 3 


1.08 x 10- 2 




11-point 


0.288 


2.115 


3.32 x 10~ 4 


4.71 x 10~ 3 




\i— int. spline 


0.415 


3.256 


4.21 x 10~ 16 


1.65 x 10~ 15 




r— int. spline 


0.238 


1.382 


1.28 x 10" 4 


1.99 x 10~ 3 


Logarithmic 


17-point 


0.665 


2.657 


1.85 x 10~ 3 


8.70 x 10" 3 




11-point 


0.345 


1.587 


5.35 x 10~ 4 


3.41 x 10" 3 




/i— int. spline 


0.528 


2.363 


1.85 x 10~ 16 


7.91 x 10- 16 




r— int. spline 


0.259 


0.999 


1.66 x 10" 4 


1.18 x 10~ 3 


Claret 


17-point 


0.190 


0.625 


8.27 x 10" 5 


5.05 x 10~ 4 




11-point 


0.175 


0.558 


4.57 x 10" 5 


2.17 x 10~ 4 




/i— int. spline 


0.184 


0.585 


4.48 x 10~ 16 


2.62 x 10~ 15 




r— int. spline 


0.170 


0.550 


1.81 x 10~ 5 


5.48 x 10~ 5 



Note. — Average and maximum values of relative rms residual a and relative flux excess 
\AF/F\ evaluated for BVRI profiles of full range of Kurucz model atmospheres (see text). 
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Table 2. Linear limb darkening coefficients measured in eclipsing binaries 



Binary 




B 


Measured u 
V 


Ft 


B 


Model uq 
V 


R 


B 


Model u r 
V 


R 


Ref. 


BS Dra 


A 


0.64 ± 0.04 


0.50 ±0.03 




0.707 


0.605 




0.651 


0.539 




1 




B 


0.64 ± 0.04 


0.50 ±0.03 




0.706 


0.605 




0.650 


0.539 






EE Peg 


A 


0.62 ±0.03 






0.660 






0.597 






2 




B 


0.75 ±0.15 






0.740 






0.689 








FS Mon 


A 


0.58 


0.58 




0.710 


0.607 




0.654 


0.540 




3 




B 


0.594 


0.591 




0.718 


0.614 




0.664 


0.548 






GG Ori 


A 


0.50 ±0.04 


0.51 ±0.03 


0.23 ±0.07 


0.594 


0.513 


0.432 


0.521 


0.446 


0.377 


4 




B 


0.50 ±0.04 


0.51 ±0.03 


0.23 ±0.07 


0.594 


0.513 


0.432 


0.521 


0.446 


0.377 




WW Cam 


A 




0.494 ±0.017 






0.601 






0.533 




5 




B 




0.499 ±0.017 






0.606 






0.538 






V459 Cas 


A 




0.487 ±0.008 






0.546 






0.479 




6 




B 




0.487 ±0.008 






0.548 






0.481 






MU Cas 


A 




0.56 ±0.07 






0.392 






0.331 




7 




B 




0.56 ±0.07 






0.399 






0.337 






WW Aur 


A 


0.616 ±0.056 


0.416 ±0.060 




0.692 


0.598 




0.625 


0.528 




8 




B 


0.512 ±0.078 


0.418 ±0.083 




0.688 


0.593 




0.618 


0.518 






RW Lac 


A 




0.55 






0.668 






0.611 




9 




B 




0.57 






0.689 






0.638 







Note. — For each eclipsing binary primary star marked by A, secondary by B. Model coefficients based on Kurucz ATLAS9 
atmospheres: uq taken from Claret (2000), u r computed in this work. Better agreeing coefficient is marked in boldface. 

References. — (1) Popper & Etzcl 1981; (2) Lacy & Popper 1984; (3) Lacy et al. 2000; (4) Torres et al. 2000; (5) Lacy ct al. 
2002; (6) Lacy, Claret, & Sabby 2004a; (7) Lacy ct al. 2004b; (8) Southworth et al. 2005; (9) Lacy et al. 2005 
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Fig. I.— Sample Kurucz ATLAS9 model intensity profile (T eff = 4250 K, log g = 4.5, 
[Fe/H] = —3.5, V band) and its different linear limb darkening fits. Original 17 points 
marked by crosses (11 points) and pluses, lines include interpolated profile / (bold solid), 
17-point fit (dotted, linear coefficient Un = 0.596), 11-point fit (short-dashed, un = 0.538), 
fit of interpolated / by \i— integration (dot-dashed, w M = 0.540), and fit of / by r— integration 
(long-dashed, u r = 0.401). Intensity shown in top panels, absolute residuals of fits from / in 
bottom panels. Left panels plotted as a function of r, right panels as a function of /x, with 
respective alternate parameters marked along the upper axes. 
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Fig. 2. — Linear limb darkening coefficient u r of solar abundance, v t = 2 km s _1 Kurucz 
model atmospheres (Kurucz 1994) plotted as a function of effective temperature for bands 
B (crosses), V (circles), R (pluses), I (triangles). 
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Fig. 3. — Relative rms residuals a and relative flux excesses AF/F for the r— integrated limb 
darkening fits from Figure 2 (left column) and for 11-point fits (right column) plotted as a 
function of effective temperature. Range of model atmospheres and symbols for photometric 
bands as in Figure 2. 
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Fig. 4. — Difference between linear limb darkening coefficients u r and tin plotted as a 
function of u n for full range of models (see text for details). 
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Fig. 5. — Difference between linear limb darkening coefficients u r and u n for full range of 
models plotted as a function of T e // separately for bands B (top left), V (top right), R 
(bottom left), / (bottom right). 
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Fig. 6. — Relative residual 7l// — 1 of linear limb darkening law I L with theoretical coeffi- 
cients u r (solid line) and Un (dashed line) compared with solar limb darkening I Q (Cox 2000) 
for four wavelengths (marked in upper right corners of the panels), plotted as a function of 
radial position r. Relative rms residuals a and relative flux excesses AF/F for both fits are 
given in each panel. 



